Christmas Bird Count scouting tool

# Create a circular buffer using the center and radius
circle <- st_buffer(st_sfc(st_point(center), crs = 4326), dist = radius_km * 1000)
bbox <- st_bbox(circle)

recent_sigtings_file <- here("Projects", "CBC_Scouting", tag, "CBC_sightings_recent.RDS")

if (file.exists(recent_sigtings_file)) {
  CBC_sightings_recent <- readRDS(recent_sigtings_file)
} else {
  # Get recent species
  sp_list <- ebirdgeo(lat = center_lat, lng = center_long, dist = radius_km, back = time_days) %>%
    pull(speciesCode)

  CBC_sightings_recent_raw <- mapply(ebirdgeo, species = sp_list, lat = center_lat, lng = center_long, dist = radius_km, back = time_days)

  # Combine in one dataframe
  CBC_sightings_recent <- bind_rows(
    lapply(CBC_sightings_recent_raw, function(df) {
      if ("exoticCategory" %in% colnames(df)) {
        # If "exoticCategory" column is present, keep it, else create it with NA values
        df <- df %>%
          mutate(exoticCategory = ifelse(!("exoticCategory" %in% colnames(df)), NA, as.character(exoticCategory)))
      } else {
        # If "exoticCategory" column is not present, create it with NA values
        df$exoticCategory <- NA_character_
      }
      return(df)
    })
  )

  saveRDS(CBC_sightings_recent, recent_sigtings_file)
}

# Map will bin observations into which week they were observed.
# First project and format date field
CBC_sightings_recent <- CBC_sightings_recent %>%
  st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
  mutate(Date = as.Date(obsDt))

# Bin dates (approximate number of weeks)
num_bins <- round(time_days / 7, 0)

# Calculate bin edges, extending the range slightly beyond the minimum and maximum dates
min_date <- min(CBC_sightings_recent$Date) - 1
max_date <- max(CBC_sightings_recent$Date)
bin_edges <- seq(min_date, max_date, length.out = num_bins + 1)

# Bin the Date variable
CBC_sightings_recent$DateBin <- cut(CBC_sightings_recent$Date, breaks = bin_edges, labels = FALSE, right = TRUE)

# Get a map tile from a map service (e.g., Stadia Alidade). This ill throw an error if API key isn't in .Renviron.
map_tile <- get_stadiamap(bbox = c(bbox[[1]], bbox[[2]], bbox[[3]], bbox[[4]]), zoom = 11, maptype = "alidade_smooth")

# Generate map
p <- ggmap(map_tile) +
  geom_sf(data = circle, fill = "transparent", color = "grey20", size = 1, inherit.aes = FALSE) +
  geom_sf(data = st_jitter(CBC_sightings_recent), size = 1.5, alpha = .9, aes(color = as.factor(DateBin)), inherit.aes = FALSE) +
  scale_color_manual(name = "Week starting:", values = rev(viridisLite::inferno(num_bins)), labels = format(bin_edges, "%Y-%m-%d")) +
  ggthemes::theme_map() +
  theme(
    legend.position = "top",
    legend.justification = "right",
    legend.background = element_rect(
      fill = "grey90",
      linewidth = 0.25, linetype = "solid",
      colour = "grey40"
    ),
    legend.key = element_rect(fill = "grey90"),
    plot.title = element_text(size = 8, face = "bold"),
    plot.subtitle = element_text(size = 6, face = "italic"),
    strip.text = element_text(size = 6)
  ) +
  facet_wrap(~comName) +
  labs(
    title = "West Hennepin CBC: recent eBird observations",
    subtitle = paste0("Data from 30 day period ending ", today)
  )

ggsave(plot = p, here("Projects", "CBC_Scouting", tag, paste0(tag, "_recent_ebird_data.pdf")), width = 8.5, height = 11)

p

recent_hotspots_sigtings_file <- here("Projects", "CBC_Scouting", tag, "CBC_hotspot_sightings_recent.RDS")
hotspots_file <- here("Projects", "CBC_Scouting", tag, "CBC_hotspots.RDS")

if (file.exists(recent_hotspots_sigtings_file) & file.exists(hotspots_file)) {
  CBC_hotspot_sightings_recent <- readRDS(recent_hotspots_sigtings_file)
  CBC_hotspots <- readRDS(hotspots_file)
} else {
  CBC_hotspots <- rebird::ebirdhotspotlist(lat = center_lat, lng = center_long, dist = radius_km) %>%
    st_as_sf(coords = c("lng", "lat"), crs = 4326)
  saveRDS(CBC_hotspots, hotspots_file)

  mapview(CBC_hotspots, label = "locName")

  CBC_hotspot_sightings_recent <- sapply(CBC_hotspots$locId, ebirdregion, simple = FALSE, back = time_days)

  # Combine in one data frame
  CBC_hotspot_sightings_recent <- bind_rows(
    lapply(CBC_hotspot_sightings_recent, function(df) {
      if ("exoticCategory" %in% colnames(df)) {
        # If "exoticCategory" column is present, keep it, else create it with NA values
        df <- df %>%
          mutate(exoticCategory = ifelse(!("exoticCategory" %in% colnames(df)), NA, as.character(exoticCategory)))
      } else {
        # If "exoticCategory" column is not present, create it with NA values
        df$exoticCategory <- NA_character_
      }
      return(df)
    })
  )

  CBC_hotspot_sightings_recent <- CBC_hotspot_sightings_recent %>%
    replace_na(list(howMany = 1))

  saveRDS(CBC_hotspot_sightings_recent, recent_hotspots_sigtings_file)
}

# Generate table with recent species observed at each hotspot (with date last documented)
species_by_hotspot <- CBC_hotspot_sightings_recent %>%
  mutate(obsDtmd = format(as.Date(obsDt), "%m-%d")) %>%
  pivot_wider(id_cols = comName, names_from = locName, values_from = obsDtmd, values_fn = max) %>%
  left_join(rebird:::tax %>% select(comName, taxonOrder)) %>%
  arrange(taxonOrder) %>%
  select(-taxonOrder)
species_by_hotspot[is.na(species_by_hotspot)] <- ""
species_by_hotspot <- species_by_hotspot %>%
  select(sort(tidyselect::peek_vars())) %>%
  select(comName, everything())

species_by_hotspot %>%
  datatable(
    options = list(pageLength = 25),
    caption = "Hotspots each species has been recently documented in"
  )

write.csv(species_by_hotspot, here("Projects", "CBC_Scouting", tag, "species_by_hotspot_last_seen.csv"))

count_hotspots_by_sp <- CBC_hotspot_sightings_recent %>%
  group_by(comName) %>%
  summarize(
    n_hotspots_observed = length(unique(locName)),
    parks_observed = paste(sort(unique(locName)), collapse = ", ")
  ) %>%
  left_join(rebird:::tax %>% select(comName, taxonOrder)) %>%
  arrange(taxonOrder) %>%
  select(-taxonOrder) %>%
  arrange(desc(n_hotspots_observed))

count_hotspots_by_sp %>%
  datatable(
    options = list(pageLength = 25),
    caption = "Number of hotspots each species has recently been observed in"
  )

write.csv(species_by_hotspot, here("Projects", "CBC_Scouting", tag, "count_hotspots_by_sp_raw.csv"))

# rare sp are those seen in fewer than 30% of hotspots
rare_sp <- count_hotspots_by_sp %>%
  filter(n_hotspots_observed <= nrow(CBC_hotspots) * .3) %>%
  pull(comName)

# unique sp are those seen in only 1 hotspot
unique_sp <- count_hotspots_by_sp %>%
  filter(n_hotspots_observed == 1) %>%
  pull(comName)

count_sp_by_hotspot <- CBC_hotspot_sightings_recent %>%
  group_by(locName) %>%
  summarize(
    n_sp = length(unique(comName)),
    n_rare_sp = length(intersect(rare_sp, unique(comName))),
    n_unique_sp = length(intersect(unique_sp, unique(comName))),
    rare_sp = paste(sort(intersect(rare_sp, unique(comName))), collapse = ", "),
    unique_sp = paste(sort(intersect(unique_sp, unique(comName))), collapse = ", ")
  ) %>%
  arrange(desc(n_sp))

write.csv(species_by_hotspot, here("Projects", "CBC_Scouting", tag, "count_sp_by_hotspot_raw.csv"))

count_sp_by_hotspot %>%
  datatable(
    options = list(pageLength = 25),
    caption = "Hotspots with total number of recently species and rare/unique species. Rare species are those seen in <=30% of hotspots."
  )